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Abstract: In order to study the mechanical behaviour 
of polar ice masses, the method of continuum me- 
chanics is used. The newly developed CAFFE model 
(Continuum-mechanical, Anisotropic Flow model, based 
on an anisotropic Flow Enhancement factor) is described, 
which comprises an anisotropic flow law as well as a fab- 
ric evolution equation. The flow law is an extension of the 
isotropic Glen's flow law, in which anisotropy enters via 
an enhancement factor that depends on the deformability 
of the polycrystal. The fabric evolution equation results 
from an orientational mass balance and includes constitu- 
tive relations for grain rotation and recrystallization. The 
CAFFE model fulfills all the fundamental principles of 
classical continuum mechanics, is sufficiently simple to 
allow numerical implementations in ice-flow models and 
contains only a limited number of free parameters. The 
applicability of the CAFFE model is demonstrated by a 
case study for the site of the EPICA (European Project 
for Ice Coring in Antarctica) ice core in Dronning Maud 
Land, East Antarctica. 

Key words: Ice, polycrystal, flow, anisotropy, continuum 
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1 Introduction 

Ice in natural land ice masses, such as polar ice sheets, 
ice caps or glaciers, consists of zillions of individual 
hexagonal crystals ("crystallites" or "grains") with a typ- 
ical diameter of millimeters to centimeters. This length 
scale stands in contrast with the size of the ice masses, 
which ranges from 100's of meters to 1000's of kilome- 
ters. It has long been known that, while the distribu- 
tion of the crystallographic axes (also known as optical 
axes or c-axes) at the surface of an ice sheet is essentially 
at random, deeper down into the ice, different types of 
anisotropic fabrics with preferred orientations of the c- 
axes tend to develop. 

Many models have been proposed to describe the 
anisotropy of polar ice. On the one end of the range 
in complexity, a simple flow enhancement factor is in- 
troduced in an ad-hoc fashion as a multiplier of the 



isotropic ice fluidity in order to account for anisotropy 
and/or impurities. This is done in most current ice-sheet 
models, often without explicitly mentioning anisotropy 
|[T9l l25l l46l . In macroscopic, phenomenological mod- 
els, an anisotropic macroscopic formulation for the flow 
law of the polycrystal is postulated. To be usable, the 
rheological parameters that enter this law must be eval- 
uated as functions of the anisotropic fabric IT4l [T5l [341 
[35l . The concept of homogenization models, also called 
micro-macro models, is to derive the polycrystalline be- 
haviour at the level of individual crystals and the fab- 
ric fl] IU E] El M E2 EB ED. As for the "high-end" 
complexity, full-field models solve the Stokes equation 
for ice flow properly by decomposing the polycrystal 
into many elements, which makes it possible to infer the 
stress and strain-rate heterogeneities at the microscopic 
scale ll27l l28l [30l [3D . A very comprehensive, up-to-date 
overview of these different types of models is given by 
Gagliardini et al. |[T3l (in this volume). However, the 
more sophisticated models are usually too complex and 
computationally time-consuming to be included readily 
in a model of macroscopic ice flow. 

Here, the Continuum-mechanical, Anisotropic Flow 
model, based on an anisotropic Flow Enhancement factor 
("CAFFE model"), shall be described. It belongs to the 
class of macroscopic models, and is laid down in detail in 
the study by Placidi et al. [42 1 (based on previous works 
by Faria flSE], Faria et al. CO, Placidi BQHm, Placidi 
and Hutter [43 1). The flow enhancement factor is taken 
as a function of a newly introduced scalar quantity called 
deformability, which is essentially a non-dimensional in- 
variant related to the shear stress acting on the basal plane 
of a crystallite, weighted by the orientation-distribution 
function which describes the anisotropic fabric of the 
polycrystal. Fabric evolution is modelled by an orienta- 
tion mass balance which accounts for grain rotation and 
recrystallization processes. 

The CAFFE model fulfills all the fundamental princi- 
ples of classical continuum mechanics (see also Placidi 
et al. [44 1), and it is a good compromise between nec- 
essary simplicity on the one hand, and consideration of 
the major effects of anisotropy on the other. In order 
to demonstrate its performance, the model is applied to 
the site of the EPICA (European Project for Ice Cor- 
ing in Antarctica) ice core in Dronning Maud Land, East 
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Quantity 


Value 










Stress exponent, n 


3 










Pre-exponential constant, Aq 


3.985 x lCT^s^Pa- 3 


(for 


T 


< 


263.15 K) 




1.916 x 10 3 s- 1 Pa" 3 


(for 


r 


> 


263.15 K) 


Activation energy, Q 


60 kJ mor 1 


(for 


r 


< 


263.15 K) 




139 kJ mor 1 


(for 


T 


> 


263.15 K) 



Table 1: Physical parameters for Glen's flow law. 



Antarctica, for which data on the ice flow as well as on 
the anisotropic fabric are available. 

2 CAFFE model 

2.1 Glen's flow law 

Let us briefly review the isotropic case, for which poly- 
crystalline ice is treated as an incompressible, viscous 
fluid. The Cauchy stress tensor T is split up according 
to 



T = —pi + S, 



1 



P 



trT, 



(1) 



where p denotes the pressure, and S is the traceless stress 
deviator (trS = 0). Due to the incompressibility, the flow 
law only determines the stress deviator S and reads 



S = 2r)D, 



(2) 



where D = sym grad v is the strain-rate tensor (symmet- 
ric part of the gradient of the velocity v), and the coeffi- 
cient ?7 is the shear viscosity (or simply the viscosity). Its 
inverse, the fluidity, can be factorized as 



1 



where 



= 2EA(T')f(a), 



itr(S2 



(3) 



(4) 



is the effective stress (square root of the second invari- 
ant of the stress deviator), and the creep function f(a) is 
given by the power law 



/(*) = ° r ' 



(5) 



(the parameter n is called "stress exponent"). Further, 
the rate factor A(T') depends on the temperature rela- 
tive to pressure melting T = T - T m + T (T: ab- 
solute temperature, T m = T — (3p: pressure melting 
point, To = 273.16 K: melting point at zero pressure, 
(3 = 9.8x 10~ 2 KMPa -1 : Clausius-Clapeyron constant) 
via the Arrhenius law 



A(T')=A e 



-Q/RT' 



(6) 



where Aq is the pre-exponential constant, Q the activa- 
tion energy and R = 8.314 Jmol -1 K _1 the universal 
gas constant. The^ow enhancement factor E is equal to 
unity for pure ice, and can be set to values deviating from 
unity in order to account roughly for effects of impurities 
and/or anisotropy (Paterson 1371). 



The isotropic flow law for ice is now obtained by in- 
serting Eq. ([3]l [with the specifications of Eqs. |5]) and 
|6]l] in the viscous flow law Q. This yields 



D = EA(T')f(a)5, 



(7) 



which is called Nye 's generalization of Glen 's flow law, 
or Glen's flow law for short (e.g., Greve and Blatter |20l . 
Hooke [23 1, Paterson OH, van der Veen |52)). Suitable 
values for the several parameters are listed in Table [T] 

2.2 Anisotropic generalization of Glen's flow law 
2.2.1 Deformation of a crystallite 

In order to derive a generalization of Glen's flow law (|7| 
which accounts for general, anisotropic fabrics of the ice 
polycrystal, we first consider the deformation of a crys- 
tallite embedded in the polycrystalline aggregate. Fol- 
lowing Placidi et al. [42], only the dominant deformation 
along the basal plane is accounted for, whereas deforma- 
tions along prismatic and pyramidal planes, which are at 
least 60 times more difficult to activate, shall be neglected 
(Fig.[l}. 



Basal 



Prismatic 



Pyramidal 



Figure 1: Basal, prismatic and pyramidal glide planes in 
the hexagonal ice crystal, sketched as a right hexagonal 
prism (Faria $9$). 

Let n be the normal unit vector of the basal plane (di- 
rection of the c-axis), then Tn is the resolved stress vec- 
tor (Fig. [2|. Note that the tensor T is interpreted as the 
macroscopic stress which does not depend on the orienta- 
tion n. It is reasonable to assume that only the stress com- 
ponent S t tangential to the basal plane (resolved shear 
stress) contributes to its shear deformation, while the 
component normal to the basal plane has no effect. 
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(Tn-n) n 



basal plane 



Figure 2: Decomposition of the stress vector into a part 
normal and a part tangential to the basal plane. 

According to Fig. [2] the decomposition of the stress 
vector reads 



Tn = (Tn • n)n + 5 t t, 



(8) 



where t denotes the tangential unit vector. Inserting the 
decomposition ([TJ of the stress tensor T readily elimi- 
nates the pressure p and leaves 



Sn = (Sn • n)n + S t t. 



(9) 



As mentioned above, deformation of the crystallite in the 
polycrystalline aggregate is attributed to the tangential 
component S t only. Since we aim at a theory which de- 
scribes the effects of anisotropy by a scalar, anisotropic 
flow enhancement factor, we define the scalar invariant 



S 2 = Sn • Sn - (Sn • n) 



(10) 



This quantity has the unit of a stress squared, and so a nat- 
ural way to non-dimensionalize it is by the square of the 
effective stress a [Eq. Q], which is also a scalar invari- 
ant. Thus, we introduce the deformability of a crystallite 
in the polycrystalline aggregate, which is loaded by the 
stress T, as 

^ (n)= 5^) =5 Sn.Sn-(Sn-^ (n) 



2 a 2 



tr(S 2 ) 



The factor 5/2 has been introduced merely for reasons of 
convenience, as it will become clear below. 



2.2.2 Flow law for polycrystalline ice 

In polycrystalline ice, the crystallites within a control 
volume (which is assumed to be large compared to the 
crystallite dimensions, but small compared to the macro- 
scopic scale of ice flow) show a certain fabric. Extreme 
cases are on the one hand the single maximum fabric, for 
which all c-axes are perfectly aligned, and on the other 
hand the isotropic fabric with a completely random distri- 
bution of the c-axes. A general fabric, which is usually in 
between these cases, can be described by the orientation 
mass density (OMD) p*(n). It is defined as the mass per 
volume and orientation, the latter being specified by the 
normal unit vector (direction of the c-axis) n G S 2 (S 2 is 



the unit sphere). Evidently, when integrated over all ori- 
entations, the OMD must yield the normal mass density 
p, which leads to the normalization condition 



J p*(n) d 2 n = p. 

S 2 



(12) 



Alternatively, an orientation distribution function (ODF) 
/*(n) can be defined as 



/» 



p*(n) 



(13) 



which is normalized to unity when integrated over all ori- 
entations. 

We use the ODF in order to define the deformability of 
polycrystalline ice by weighting the deformability of the 



crystallite (Hi, 



A = J A*{n)f*(n)d 2 n 

S 2 

s 2 

_ f Sn • Sn - (Sn • n) 2 2 
= 5 J ir^ f (n)d ^ (U) 



Note that, for the isotropic case, the ODF is /*(n) = 
l/(47r), and from Eq. (fl4| we obtain a deformability of 
A = 1 (Placidi et al. [42]). For that reason, the factor 5/2 



has been introduced in Eq. (Hi. 



UC/SM 



I 



SS/SM 



t 



t 



T 



Figure 3: Uniaxial compression on single maximum 
(UC/SM) and simple shear on single maximum (SS/SM) 
for a small sample of polycrystalline ice. Stresses are in- 
dicated as black arrows, and the single maximum fabric 
is marked by the dark-grey arrows within the ice sample. 

The envisaged flow law for anisotropic polar ice can 
now be formulated. Essentially, we keep the form of 
the Glen's flow law ([7), but with a scalar, anisotropic en- 
hancement factor E(A) instead of the parameter E, 

D = E(A)A(T')f(a)S. (15) 

The function E(A) is supposed to be strictly increasing 
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with the deformability A, and has the fixed points 



2.2.3 Inversion of the flow law 



(16) 



E{0) = E min 

(uniaxial compression on single maximum), 
E(l) = 1 

(arbitrary stress on isotropic fabric), 

= £max 

(simple shear on single maximum). 



The "hard" case ([T6|i and the "soft" case ( |T6] >3 are illus- 
trated in Fig. [3] Note also that the deformability cannot 
take values larger than .4 = 5/2 (Placidi et al. [42]). 

For the detailed form E(A) of the anisotropic en- 
hancement factor, in addition to Eq. (16 1, we demand 
that the function is continuously differentiable, that is, 



E e C 1 [0, §]. Moreover, Azuma |2| and Miyamoto 
have experimentally verified that the enhancement factor 
depends on the Schmid factor (resolved shear stress) to 
the fourth power, that is, on the square of the deformabil- 
ity A. This yields 



E n 



(1 - E^)A l 



21 l-E„ 



E(A) = < 



<A< 1, 



(17) 



4^ 2 (^ ma x-l) + 25-4£ n 



21 



1< A< 



(for details see Placidi et al. 11421 ). Several investigations 
(e.g. Budd and Jacka [3|, Pimienta et al. 1391 . Russell- 
Head and Budd [45 1) indicate that the parameter E max 
(maximum softening) is approximately equal to ten. The 
parameter E m [ D (maximum hardening) can be realisti- 
cally chosen between zero and one tenth, a non-zero 
value serving mainly the purpose of avoiding numerical 
problems. The function ( 17 1 is shown in Fig. [4] 
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Figure 4: Anisotropic enhancement factor E(A) as a 
function of the deformability A according to Eq. (|77|), for 



As long as the creep function f(a) is given by the power 
law the anisotropic flow law ( p~5] > can be inverted an- 
alytically. We find 

S = [E(AT 1/n [A(T')]- 1/n d- {1 - 1/n) D, (18) 

where 

(19) 



d = V|tr(D2) 



is the effective strain rate. The deformability A also 
needs to be expressed by strain rates instead of stresses 



[see Eq. (14i]. In analogy to Eq. (|9J, we consider the 
resolved strain-rate vector Dn in a crystallite in the poly- 
crystalline aggregate, and decompose it according to 



Dn= (Dnn)n + At, 



(20) 



where D t is the resolved shear rate in the basal plane (see 
also Fig. [2]). As in Eq. ( 10 1, we define the scalar invariant 



D 2 t = Dn • Dn - (Dn • n) 2 



(21) 



Owing to the collinearity of the tensors S and D [see 
Eqs. ( 15 1 and ([l8j], the deformability of a crystallite in 
the polycrystalline aggregate [Eq. < |TT| >] can be readily ex- 
pressed by D t and d, 

s 5 fl t 2 (n) Dn- Dn - (Dn ■ n) 2 

^ = 2-^ = 5 — *m ■ (22) 

and the deformability of polycrystalline ice [Eq. ( fl4} ] 
yields 



A = J A k {n)f{n)A 2 n 

s 2 



« / D - D °- D ^-» 2 r( n)dV ,23) 



This completes the inversion of the anisotropic flow law. 



2.3 Evolution of anisotropy 
2.3.1 Orientation mass balance 



E„ 



10 and E,, 



0. 



The anisotropic flow law in the form ( 15 1 or ( 18 i needs 
to be complemented by an evolution equation for the 
anisotropic fabric. This is done by formulating the ori- 
entation mass balance for the OMD /9*(n). 
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Figure 5: Orientation transition rate u*(n) on the unit 
sphere S 2 . 

We are not going to enter the detailed formalism of 
orientation balance equations here (see Placidi et al. Il42l . 
and references therein). Instead, we rather motivate the 
form of the orientation mass balance by generalizing the 
ordinary mass balance. The difference is that, in addi- 
tion to the dependencies on the position vector x and the 
time t, the density and velocity fields also depend on the 
orientation vector n £ S 2 , which is indicated by the no- 
tation p*(n) and v*(n). The velocity, which describes 
motions in the physical space, is complemented by an ori- 
entation transition rate u*(n), which describes motions 
on the unit sphere, that is, changes of the orientation due 
to grain rotation (Fig.[5|l. Also, an orientation flux q*(n) 
is considered, which allows redistributions of the OMD 
due to rotation recrystallization (polygonization). Conse- 
quently, the orientation mass balance reads 

dp* 

+ div (p*v*) + div s2 (p*u* + q*) = p*T* . (24) 

at 

The first two terms on the left-hand side are straightfor- 
ward generalizations of the terms in the ordinary mass 
balance. The third term on the left-hand side is the equiv- 
alent of the second term for the orientation transition rate 
u*(n) and the orientation flux q*(n), where div ,52 is the 
divergence operator on the unit sphere. On the right-hand 
side, a source term appears which allows that certain ori- 
entations can be produced at the expense of others. The 
quantity r*(n) is therefore called the orientation produc- 
tion rate. Physically, it describes migration recrystalliza- 
tion and all other processes in which the transport of mass 
from one grain, having a certain orientation, to another 
grain, having a different orientation, cannot be neglected. 

In the following, we will make the reasonable assump- 
tion that the spatial velocity does not depend on the ori- 
entation, that is, v*(n) = v. Therefore, the orientation 
mass balance p4| ) simplifies to 

dp* 



dt 



+ div (p*v) + div S 2 (p*u* + q*) = p*T*. (25) 



Integration over S 2 (all orientations) gives the classical 
mass balance 



with the use of the Gauss theorem and the mass- 
conservation requirement 

J p*(n)r*(n)d 2 n = 0. (27) 

s 2 

In order to solve the orientation mass balance (25) , 
constitutive relations for the orientation transition rate 
u*(n), the orientation flux q*(n) and the orientation pro- 
duction rate T* (n) need to be provided as closure condi- 
tions. 



2.3.2 Constitutive relation for the orientation transi- 
tion rate 

As mentioned above, the orientation transition rate corre- 
sponds physically to grain rotation. Since grain rotation 
is induced by shear deformation in the basal plane, we 
argue that it is controlled by the resolved shear rate D t t 
[Eq. <|20]>], and use the relation 

u*(n) = -tAt+Wn 

= i- [(Dn-n)n- Dn] + Wn (28) 

(see, e.g., Dafalias 0). The parameter l is assumed to 
be a positive constant. The additional term Wn with the 
spin tensor W = skw grad v (skew-symmetric part of the 
gradient of the velocity v) describes the contribution of 
local rigid-body rotations. 

In the special case 1 = 1, the basal planes are material 
area elements, that is, they carry out an affine rotation. 
However, due to geometric incompatibilities of the de- 
formation of individual crystallites in the polycrystalline 
aggregate, an affine rotation is not plausible, and we ex- 
pect realistic values of u to be less than unity. 

2.3.3 Constitutive relation for the orientation flux 

The orientation flux is supposed to describe rotation re- 
crystallization (polygonization). Following the argumen- 
tation by Godert ifTBI . it is modelled as a diffusive pro- 
cess, 

q*(n) = -Agrad s2 [/(n)^*(n)], (29) 

where the parameter A > is the orientation diffusivity, 
gradg2 is the gradient operator on the unit sphere, and the 
"hardness" H* (n) is a monotonically decreasing function 
of the crystallite deformability A*(n) [see Eq. (111]. A 
simple choice for the hardness function would therefore 
be 



H*(n) 



1 



(30) 



A*(n) + e : 

the offset e < 1 being introduced in order to prevent a 
singularity for A* — 0. However, recent results by Du- 
rand et al. suggest that rotation recrystallization is an 
isotropic process not affected by the orientation. In this 
case, the choice 

H*(n) = 1 (31) 



dp 
dt 



div (pv) = 0, 



(26) 



is indicated, which renders Eq. ( 29 1 equivalent to Fick's 
laws of diffusion on the unit sphere. 
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2.3.4 Constitutive relation for the orientation pro- 
duction rate 

The driving force for the orientation production rate, 
which models essentially migration recrystallization, are 
macroscopic deformations of the polycrystal, which can 
be more easily followed on the microscopic scale by 
grains oriented favourably for the given deformation. 
Therefore, it is reasonable to assume that the orientation 
production rate for a certain orientation n is related to the 
crystallite deformability _4*(n) [Eqs. 
CAFFE model, the linear relation 



Hi, (22 1]. In the 



r*(n) =T[A*(n)-A] 



(32) 



is proposed. Subtraction of the polycrystal deformabil- 
ity A is required in order to fulfill the mass-conservation 
condition §2J\ . The parameter T is assumed to be pos- 
itive, which guarantees a positive mass production for 
favourably oriented grains, and a negative production for 
unfavourably oriented grains (Pig. [6j. 




an abrupt change towards a single maximum, which pre- 
vails below 2040 m depth. Tendencies of secondary or 
multiple maxima are observed at several depths. The 
complete data set and a detailed interpretation will be pre- 
sented elsewhere (Hamann et al. J22|). 

The location of the EDML site on a flank (rather than 
a dome like most other ice cores) allows deriving a one- 
dimensional flow model based on the shallow-ice approx- 
imation (Hutter [24|, Morland [33]), with which the per- 
formance of the CAFFE model can be tested. We define a 
local Cartesian coordinate system such that Kohnen Sta- 
tion is located at the origin, the ir-axis points in the 260° 
(WSW) direction, the y-axis in the 170° (SSE) direction, 
and z (depth) points vertically downward (Fig. [7}. 



-74°48'-r 



-74-51'- 




x (260°) 

-75°03' 



-75°06' 



0°00' * y (170°) 0°30' 
5 10 km 



Figure 6: 
Eq. 02). 



Orientation production rate according to 



The CAFFE model is now formulated completely. 



Equation (15i is the actual flow law, which replaces its 
isotropic counterpart (j7]). Anisotropy enters via the en- 
hancement factor E(A) [Eq. \l7\], which depends on the 
deformability A defined in Eq. (|14|i. Computation of the 
deformability requires knowledge of the orientation mass 
density p*, which is governed by the evolution equation 
and the constitutive relations (EHl, d29l) and ((32l). 



3 Application to the EDML ice core 

3.1 Methods 

We have developed a one-dimensional flow model, in- 
cluding the CAFFE model, for the site of the EPICA ice 
core at Kohnen Station in Dronning Maud Land, East 
Antarctica ("EDML core", 75°00'06"S, 00°04'04"E, 
2892 meters above sea level; see http://www.awi- 
bremerhaven.de/Polar/Kohnen/). For this core with an 
overall length of 2774 m, preliminary fabric data are 
available from 50 m until 2570 m depth (I. Hamann, pers. 
comm. 2006). The fabric is essentially isotropic down to 
approximately 600 m depth, and shows a gradual tran- 
sition to a broad girdle fabric between 600 and 1000 m 
depth. Further down, the girdle fabric narrows until ap- 
proximately 2000 m depth. The fabric then experiences 



Figure 7: Local coordinate system for the EDML site. 
Underlaid topography map by Wesche et al. H53V . 

According to the topographical data by Wesche et al. 
Il53l . the x-axis is approximately aligned with the down- 
hill direction, and the gradient of the free surface eleva- 
tion, h, is 



dh 
dx 



-9 x 10~ 4 ± 10%, 



dh 
dy 



0. 



(33) 



Thus, in the shallow-ice approximation, the only non- 
vanishing bed-parallel shear-stress component is T xz 
(= S xz ), given by 



T xz = S xz = pgz 



dh 

dx ' 



(34) 



where g is the acceleration due to gravity. Combination 
with the x-z-component of the Glen's flow law (|7]) yields 
the isotropic horizontal velocity, 



dh 

dx 



H 



J AfT')^ 1 zdz (35) 



(e.g., Greve [181, Greve and Blatter [20|), where H is the 
ice thickness, the rate factor A(T') and stress exponent 
n are chosen as listed in Sect. |2.1| and the enhancement 
factor E has been set to unity. Similarly, for anisotropic 
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Vertical strain-rate (a -1 x 10~ 5 ) 
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Figure 8: Dansgaard-Johnsen type distributions of the vertical strain rate (left panel) and the temperature at the EDML 
site (right panel). The depth of the kinks is at two-thirds of the local ice thickness. The strain rate at the surface has been 
chosen such that the downward vertical velocity equals the accumulation rate, and the surface and basal temperatures 
match the ice-core data. 



conditions and the corresponding flow law (15 1, the hori- 
zontal velocity is 



dh 



H 



v, : = -2pffg; / E(A) A(T) a"" 1 zdz, (36) 



with the enhancement factor function E(A) of Eq. ( 17 1. 
Note that no-slip conditions have been assumed at the ice 
base, that is, v x (z — H) = 0. 

The unknowns in Eq. ([36) are the normal deviatoric 



stresses (S xx , S yy , S zz ) which are required together with 



the shallow-ice shear stress (34i for computing the de 



formability A by Eq. ( 14 1, and then the enhancement fac- 
tor E(A) by Eq. ( 17 i. The normal deviatoric stresses are 



computed by application of the inverse anisotropic flow 
law ( p"8"} with the deformability in the form ((23). The 
latter is evaluated with the calculated shallow-ice defor- 
mations and an assumed vertical strain rate D zz in the 
form of a Dansgaard-Johnsen type distribution |7], which 
consists of a constant value of D zz from the free surface 
down to two thirds of the ice thickness, and a linearly 
decreasing value of D zz below. A similar distribution is 
employed for the temperature profile (see Fig. [8). We also 
assume extension in the a; -direction only, so that the only 
non-zero horizontal strain rate entering the evaluation of 
Eq. ( f23) is D xx — —D zz . The vertical velocity v z results 
from integrating the prescribed vertical strain rate D zz , 
which gives a linear/quadratic profile (e.g., Greve et al. 

ED). 

For the ODF, we use the preliminary data of the 
EDML fabric described above. However, since during 
the drilling process the orientation of the core is not fixed, 
the horizontal orientation of the non-circularly symmet- 
ric girdle fabric (between approximately 600 and 2040 m 



depth) relative to our coordinate system, i.e., the direc- 
tion of ice flow, is unknown. For this reason, we need to 
assign an orientation for the fabric when computing the 
enhancement factor. We consider two limiting cases by 
rotating the initial data such that the girdle fabric at all 
depths is aligned with the x-axis (case "R13") and with 
the y-axis (case "R23"), respectively. This is illustrated 
in Fig. [9] 




Figure 9: Sketch of the rotation of the girdle fabrics in 
order to align with the x-axis (case "R13") and with the 
y-axis (case "R23") in the Schmidt projection. 

At the surface, we assume isotropic conditions, so that 
As = 1, and for the maximum softening and hardening 
parameters, we use the values -E max = 10 and E nlin = 0, 
respectively. 

In a second step, we attempt at solving the fabric evo- 



lution equation ( 25 1. For the lack of better knowledge, we 
neglect recrystallization, that is, we set A = and T = 
in the constitutive relations ((29) and d32), respectively. 



By allowing only a dependency of the orientation mass 
density p* on the vertical coordinate z (one-dimensional 
steady-state problem) and on the orientation n, the orien- 
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tation mass balance ( |25| l yields an equation which gov 
erns the fabric evolution along the EDML ice core, 

dp* 



dz 



v z + di(p*u*)=0, 



(37) 



where the orientational gradient operator di and the ori- 
entation transition rate u*, respectively, read in index no- 
tation as 

di = UiUj ^— (38) 



drii 

and, due to Eq. ( |28"j ), 

u* = iD hk n h n k ni 



dn,; 



W l3 n 3 . (39) 



With Eq. ( pS) , and by inserting the constitutive relation 
37) , it follows 

v z + u*d t p* + p*diU* 



Eq. (|39) in Eq. ((37), it follows 

dp* 
dz 
dp* 



dz 



+ (Wij - iDij)njdip* 
+ Zip* D hk n k n h = 0. 



(40) 



We assume that the local flow field consists of vertical 
compression with the compression rate (negative vertical 
strain rate) e = —dv z /dz according to Fig. [8] horizontal 
extension in a; -direction, and the horizontal, bed-parallel 
shear rate 

7 = fe = 2w i^- 4)A(T ' )an " z (41) 

that results from Eq. ( |36] >. The velocity gradient L = 
grad v then reads 




7 




(42) 



Consequently, we obtain for the strain-rate tensor D and 
the spin tensor W 



and 



D 



W 
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(43) 



(44) 



With these expressions and the introduction of spherical 



coordinates, Eq. (40 1 reduces to 



.dp* 
dz 



3tp* [e(2 sin 2 9 cos 2<p - 
+ 27 sin 28 cos <p\ 

dp* 



1 -3 cos 20) 



d<p 
dp* 

~d~0 



eisin2iy9 + 7(— 1 + l) 



tan 9 



-ie sin 29 (cos 2ip + 3) 



+ 7(1 — 1 cos 29) cos <p 



0, 



(45) 



where 9 and ip are the polar angle (co-latitude) and az- 
imuth angle (longitude), respectively. Note that, due to 
Eq. ( pT) , the shear rate 7 depends on the fabric via the 
deformability A. 

The shear flow at the EDML station leads to the 
transport of ice particles over significant horizontal dis- 
tances. Huybrechts et al. [25 1 estimate, based on three- 
dimensional flow modelling, that particles at 89% depth 
of the core originate from ps 184 km upstream. This is 
not taken into account in our spatially one-dimensional 
model. However, the variation of the shear upstream of 
the drill site is likely small due to the small variation of 
the surface gradient (Fig. (7), so that the error resulting 
from the neglected horizontal inhomogeneity should be 
limited. 



In this study, we restrict the solution of Eq. ((45) to 
the simplified case of a transversely isotropic (circularly 
symmetric) fabric, so that the OMD p* is only a func- 
tion of the depth z and the polar angle 9. Then Eq. (45 1 
becomes, after integration over the azimuth angle <p, 



^ z -%,^28 
dz d8 



3ip*e(l + 3cos2(9) = 0. 



(46) 



Equation (46i is solved by using a finite-difference dis- 



cretization with the parameter 1 = 0.6. 



3.2 Results 



Figure 10 shows the variation of the enhancement fac- 
tor, the ice fluidity and the horizontal velocity along the 
ice core, computed with the ODF based on the fabric 
data described above. For both limiting cases R13 and 
R23, the enhancement factor is close to unity in the up- 
per 600 m, which reflects the nearly isotropic fabrics in 
that part of the EDML core. Further down, in the girdle 
fabric regime, the case R13 is characterized by a moder- 
ate increase of the enhancement factor to an average value 
of about two, whereas the case R23 exhibits a strong de- 
crease of the enhancement factor to values close to zero. 
This demonstrates clearly that the girdle fabrics produce 
a significantly different mechanical response depending 
on the orientation relative to the ice flow. Case R23 is 
probably closer to reality, because in the girdle fabric 
regime above 2000 m depth the deformation is essentially 
pure shear (vertical compression, horizontal extension in 
x-direction only). For this situation, a simple "deck-of- 
cards" model illustrates that the c-axes turn away from 
the £-axis and towards the z-axis, whereas nothing hap- 
pens in ^-direction, so that in the Schmidt projection a 
concentration perpendicular to the a>axis (flow direction) 
results. 

Below 2000 m depth, where the fabric switches to a 
single maximum, the difference between the cases R13 
and R23 essentially vanishes. The crystallite basal planes 
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Figure 10: Variation of the enhancement factor (left panel), the ice fluidity (middle panel) and the horizontal velocity 
(right panel) along the EDML ice core. "Data R13" and "Data R23" represent the solutions obtained with the measured 
girdle fabrics rotated to align with the x- and y-direction, respectively, and "Isotropy" represents isotropic conditions. 



are favourably oriented for the now prevailing simple- 
shear deformation, which leads to large deformabilities. 
Consequently, the enhancement factor shows a sharp in- 
crease to a maximum value of about eight, which is close 
to the theoretical maximum of E max = 10. 

The variabilities of the enhancement factor and the ef- 
fective stress, as well as the increase of the temperature 
with depth, contribute to the fluidity profiles shown in 
Fig. [Top . Since the fluidity is very small above 2000 m 
depth and increases only further down, the difference be- 
tween the cases R13 and R23 in absolute values is sur- 
prisingly small. At 2563 m depth, the fluidity is about 
200 times higher than the fluidity at 1000 m depth for the 
case R23 due to the counteracting contributions from the 
favourably oriented c-axes, the higher temperature and 
the smaller effective stress. The latter is somewhat sur- 
prising; it is caused by the normal deviatoric stresses S xx 
and S zz , which decrease strongly below 2000 m depth 
and outweigh the influence of the increasing shear stress 
S xz in the effective stress. 

Owing to the large enhancement factors close to the 
bottom, the anisotropic flow law predicts significantly 
larger horizontal velocities compared to the isotropic flow 



Let us now turn to the simulation in which the fabric 



law for the entire depth of the ice core (Fig. 10 ;). At the 
surface, the anisotropic horizontal velocities are by ap- 
proximately a factor 3.5 larger than their isotropic coun- 
terparts, and the absolute value of sa 0.7 ma -1 agrees 
very well with measurements (H. Oerter, pers. comm. 
2005; Wesche et al. [53 1). The difference between the 
cases R13 and R23 amounts to rj 10%, the larger val- 
ues being obtained for the case R23 owing to the slightly 
larger enhancement factors below 2000 m depth. Inter- 
estingly, these differences show that the fabrics are not 
perfectly transversely isotropic below 2000 m depth, even 
though they are very close to the single-maximum type. 



evolution is computed by solving Eq. (46i for a trans- 
versely isotropic fabric. Although this assumption is not 
consistent with the observed girdle fabric between ap- 
proximately 600 and 2000 m depth and is therefore a 
gross simplification, it is interesting to study the mechan- 
ical response of such a simplified system and the differ- 
ences to the ice flow resulting from applying the mea- 
sured fabrics. 

Figure [TTfr shows the comparison between the en- 
hancement factors resulting from the computed, trans- 
versely isotropic fabric (which will be referred to as 
"modelled enhancement factor" in the following) and 
from the fabric data. Evidently, the agreement is good 
despite the assumption of transverse isotropy. Down to 
1800 m depth, the modelled enhancement factor lies in 
between the cases R13 and R23, which are the limiting 
cases for the orientation of the measured girdle fabric 
with respect to the ice-flow direction. Between 1800 and 
1900 m depth, the modelled enhancement factor is very 
close to the low values of the case R23, for which the 
girdle fabric is aligned perpendicular to the flow direc- 
tion. Below 2000 m depth, the sharp increase is also well 
reproduced; however, the maximum of the modelled en- 
hancement factor is more pronounced and lies closer to 
the bottom than for the cases R13 and R23. 

For that reason, the modelled enhancement factor leads 
to larger near-basal shear rates than the enhancement fac- 
tor based on the cases R13 and R23. Consequently, the 
horizontal velocity resulting from the modelled enhance- 
ment factor is larger by about a factor two than the ve- 
locities for the cases R13 and R23 (Fig. [Tip) . At the sur- 
face, a value of ~ 1.5 ma^ 1 is reached, which is twice 
the measured surface velocity. This highlights the great 
sensitivity of the ice dynamics to the processes near the 
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Figure 11: Variation of the enhancement factor (left panel) and the horizontal velocity (right panel) along the EDML ice 
core. 

R13", "Data R23" and "Isotropy" see the caption of Fig. [TO 



"Model" represents the solutions based on the fabric evolution equation {46) for transverse isotropy. For "Data 



bottom, which are most difficult to model precisely. Be- 
side the assumption of transverse isotropy, a weak point 
in that context is the neglection of recrystallization pro- 
cesses, which are expected to become important for the 
fabric evolution in the lower part of the ice core. This 
point requires further attention. 

4 Conclusions 

The newly developed CAFFE model (Continuum- 
mechanical, Anisotropic Flow model, based on an 
anisotropic Flow Enhancement factor), which comprises 
an anisotropic flow law as well as a fabric evolution equa- 
tion, was presented in this study. It is a good compro- 
mise between physical adequateness and simplicity, and 
is therefore well suited for being used in flow models of 
ice sheets and glaciers. 

The CAFFE model was successfully applied to the site 
of the EDML ice core in East Antarctica. Two different 
methods were employed, (i) computing the anisotropic 
enhancement factor and the horizontal flow based on fab- 
rics data, and (ii) solving the fabric evolution equation 
under the simplifying assumption of transverse isotropy. 
Method (i) demonstrated clearly the importance of the 
anisotropic fabric in the ice column for the flow veloc- 
ity, and better agreement with the measured surface ve- 
locity was achieved compared to an isotropic computa- 
tion. The anisotropic enhancement factor produced with 
method (ii) agreed reasonably well with that of method 
(i), despite the fact that the measured fabric is not trans- 
versely isotropic in large parts of the ice core. 

A solution of the fabric evolution equation ( |4"5j ) for 
the EDML ice core without the assumption of trans- 
verse isotropy has been presented elsewhere (Seddik et al. 
11481 ). Further, the CAFFE model has already been im- 



plemented in the three-dimensional, full-Stokes ice-flow 
model Elmer/Ice (Seddik |47), Seddik et al. |g9)) in or- 
der to simulate the ice flow in the vicinity within 100 km 
around the Dome Fuji drill site (Motoyama 061 ) in cen- 
tral East Antarctica. 
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